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We explore the mechanism responsible for the ergodicity breaking in systems with long-range 
forces. In thermodynamic limit such systems do not evolve to the Boltzmann-Gibbs equilibrium, 
but become trapped in an out-of-equilibrium quasi-stationary-state. Nevertheless, we show that if 
the initial distribution satisfies a specific constraint — a generalized virial condition — the quasi- 
stationary-state is very close to ergodic and can be described by Lynden-Bell statistics. On the 
other hand if the generalized virial condition is violated, parametric resonances are excited, leading 
to chaos and ergodicity breaking. 
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Statistical mechanics of systems in which particles in- 
teract through long-ranged potentials is fundamentally 
different from the statistical mechanics of systems with 
short-range forces pQ. In the latter case, starting from 
an arbitrary initial condition (microcanonical ensemble) 
systems evolve to a thermodynamic equilibrium in which 
particle distribution functions are given by the usual 
Boltzmann-Gibbs statistical mechanics [3J. The state 
of thermodynamic equilibrium does not depend on the 
specifics of the initial distribution, but only on the global 
conserved quantities such as energy, momentum, angu- 
lar momentum, etc. The situation is very different for 
systems in which particles interact through long-range 
potentials, such as gravity or unscreened Coulomb in- 
teractions [3HS]. In this case, it has been observed 
in numerous simulations that these systems do not re- 
lax to thermodynamic equilibrium, but become trapped 
in a quasi-stationary state (qSS), the lifetime of which 
diverges with the number of particles [3J [BHH]- The 
distribution functions in this quasi-stationary state do 
not obey the Boltzmann-Gibbs statistical mechanics — 
and in particular, particle velocities do not follow the 
Maxwell-Boltzmann distribution, but depend explicitly 
on the initial condition. It has been an outstanding chal- 
lenge of statistical mechanics to quantitatively predict 
the final stationary state reached by systems with un- 
screened long-range forces, without having to explicitly 
solve the iV-body dynamics or the collisionless Boltz- 
mann (Vlasov) equation. 

Some 40 years ago Lynden-Bell (LB) proposed a gener- 
alization of the Boltzmann-Gibbs statistical mechanics to 
treat systems with long-range interactions [9 . Lynden- 
Bell's construction was based on the Boltzmann count- 
ing, but instead of using particles, LB worked directly 
with the levels of the distribution function. The motiva- 
tion for this approach was the observation that dynami- 
cal evolution of the distribution function for systems with 
long-range interactions is governed by the Vlasov equa- 
tion |10j . This equation has an infinite number of con- 



served quantities, Casimirs — any local functional of the 
distribution function is a Casimir invariant of the Vlasov 
dynamics. In particular if the initial distribution func- 
tion is discretized into levels, the volume of each level 
must be preserved by the Vlasov flow. For an initially 
one-level distribution function, Vlasov dynamics requires 
that the phase space density does not exceed that of the 
initial distribution — one-particle distribution function 
over the reduced phase space (/x-space) evolves as an in- 
compressible fluid. Using this constraint in a combina- 
tion with the Boltzmann counting, LB was able to de- 
rive a coarse-grained entropy, the maximum of which he 
argued should correspond to the most-probable distri- 
bution — the one that should describe the equilibrium 
state. Numerous simulations, however, showed that, in 
general, Lynden-Bell statistics was not able to account 
for the particle distribution in self-gravitating systems, 
and the theory has been abandoned in the astrophysi- 
cal context. Recently, however, Lyndon-Bell's work has 
been rediscovered by the Statistical Mechanics commu- 
nity, which showed that for some systems, specifically the 
widely studied Hamiltonian Mean Field Model (HMF), 
Lynden-Bell's approach could make reasonable predic- 
tions about the structure of the phase diagram |llj . The 
fundamental question that needs to be addressed is: Un- 
der what conditions can Lynden-Bell statistics be used 
to accurately describe systems with long-range interac- 
tions? This will be the topic of the present Letter. 

To be specific, we will study the HMF model [T|, which 
has become a test bench for theories of systems with long- 
range forces. However, our results and methods are com- 
pletely general and can be applied to other systems, such 
as self-gravitating clusters or confined non-neutral plas- 
mas. The HMF model consists of N particles restricted 
to move on a circle of radius one. The dynamics is gov- 
erned by the Hamiltonian 
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where the angle 9i is the position of i'th particle and 
Pi is its conjugate momentum |llH13j . The macroscopic 
behavior of the system is characterized by the magne- 
tization vector M = (M x ,M y ), where M x = (cos 9), 
My = (sm0), and (• • •) stands for the average over all 
particles. The Hamilton's equations of motion for each 
particle reduce to 



= -M x {t) Bin 0i(t) +M y (t) coB0i(t). 
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Since the Hamiltonian does not have explicit time depen- 
dence, the average energy per particle, 
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is conserved. 

The failure of LB theory in the astrophysical context 
was attributed to incomplete relaxation, lack of good 
mixing, or broken ergodicity |14j . The mechanisms be- 
hind this failure have not been elucidated. On the other 
hand, it has been recently observed that if the initial 
distribution is virialized — satisfies the virial condition 
- LB's approach was able to quite accurately predict 
the stationary state of gravitational and Coulomb sys- 
tems .'! (i . Unfortunately, the virial theorem can be 
derived only for potentials which are homogeneous func- 
tions. This is not the case for the HMF model. Neverthe- 
less, the fact that LB theory seems to apply under some 
conditions makes one wonder if such conditions can be 
found for arbitrary long-range potentials, which are not 
in general homogeneous functions. 

To answer the questions posed above, we note that if 
the initial distribution is virialized, macroscopic oscilla- 
tions of observables should be diminished. On the other 
hand, if the system is far from virial, the mean-field po- 
tential that each particle feels will undergo strong oscil- 
lations. It is then possible for some particles to enter in 
resonance with the oscillations of the mean-field, gain- 
ing large amounts of energy. The parametric resonances 
will result in the occupation of regions of the phase-space 
which are highly improbable, from the point of view of 
Boltzmann-Gibbs or LB statistics [IB]. Furthermore, res- 
onant particles will take away energy from collective os- 
cillations producing a form of non-linear Landau damp- 
ing [15]. After some time, macroscopic oscillations will 
die out and each particle will feel only the static mean- 
field potential. From that point on, particle dynamics 
will become completely regular, with no energy exchange 
possible between the different particles. The particles 
which have gained a lot of energy from the parametric 
resonances will be trapped forever in the highly improb- 
able regions of the phase-space, unable to thermalize with 
the rest of the system. 

To see how the theoretical picture advocated above 
can be applied to the HMF, we first derive a generalized 
virial condition for this model. For simplicity we will 



consider initial distributions of the "water-bag" form in 
(0,p). Without loss of generality, we choose a frame of 
reference where (9) = and (p) = 0. The one-particle 
initial distribution function then reads 
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where O is the Heaviside step function, and \9q\ and |po| 
are the maximum values of angle and momentum, re- 
spectively. Note that from symmetry, M y (t) = at all 
times. When the dynamics starts, the mean-squared par- 
ticle position will evolve with time. We define the enve- 
lope of the particle distribution as 9 e (t) = y^3(9 2 ), so 
that 9 e (0) = 0q. We next differentiate 9 e (t) twice with 
respect to time to obtain the envelope equation of mo- 
tion, 
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Using the conservation of energy, (p 2 ) = 2u + M 2 (t) — 1. 
To calculate, (99), we use the equation of motion for 9. 
Supposing that the distribution of angles remains close 
to uniform on the interval [— 9 e , 9 e ], we obtain 
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= M x (t) cos 9 e (t)~ M 2 (t), 
where the magnetization M x (t) is 
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Neglecting the correlations between positions and veloci- 
ties, (9p) = 0, we finally obtain a dynamical equation for 
the envelope 



(2u + M x (t) cos9 e (t) - 1), 



(8) 



where u = p 2 /Q + (1 - M 2 )/2 and M = sm(0 o )/0 o . The 
generalized virial condition is defined by the stationary 
envelope, 9 e = 0, which means that along the curve 



(2m - l)6»o + sin O cos 9 = 0. 
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magnetization remains approximately invariant. In Fig. 
1 we plot Eq.Q in the Mq — u plane and compare it 
with the full molecular dynamics simulation of the HMF 
model. As can be seen, agreement between the theory 
and simulation is excellent. 

Along the generalized virial condition curve, Eq. |H, 
the magnetization — and, therefore, the mean-field po- 
tential acting on each particle — of the HMF model has 
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FIG. 1. Phase diagram of the HMF model obtained using 
the molecular dynamics simulations. Solid curve is the line 
of first order phase transitions separating paramagnetic and 
ferromagnetic phases. This line extends up to Mo = 0.6, after 
which point the order of the phase transition, shaded region 
(green line), becomes unclear, with strong dependence on the 
initial conditions and various reentrant transitions occurring 
in this region. Dashed curve is the generalized virial condi- 
tion, Eq.|9]|. Along this curve oscillations of the envelope are 
suppressed. Diamonds are the results of simulation. Start- 
ing with the initial energy and magnetization along the virial 
curve, diamonds show the final magnetization to which the 
system relaxes. For points along this curve, the final magne- 
tization is almost identical to the initial one. Note that the 
generalized virial curve terminates at Mo = 0.34 slightly be- 
low the phase transition line. This small difference, however, 
is sufficient to invalidate the Lynden-Bell theory, which for 
Mo — 0.4 predicts a second order phase transition, while the 
simulations show that the phase transition is of first order 
|17| . Points (A), (B), and (C) correspond to the initial con- 
ditions for the distribution functions shown in Fig. 2. The 
Poincare sections of the test particle dynamics for the initial 
conditions described by the points (B) and (T) are shown in 
Fig. 3. Finally, we note that since the stationary distribution 
must satisfy the virial condition and the energy is conserved, 
Eq.Q allows us to predict the magnetization to which the 
system will evolve for initial conditions lying inside the ferro- 
magnetic region, see the arrows for points (A) and (C). 



only microscopic oscillations and the parametric reso- 
nances are suppressed. Under these conditions, we ex- 
pect that LB theory will be valid. The coarse grained 
entropy within the LB approach is given by 
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where 770 — 1/ 46*oPo [IB] • Maximizing this entropy under 
the constraints of energy and particle conservation, we 
obtain the equilibrium distribution function 
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The Lagrange multipliers fj, and /3 are determined by par- 
ticle and energy conservation, 
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respectively, and the magnetization by the self- 
consistency requirement, 



dp d6 cos 0f(p,0) = M x . 



(14) 



Solving these equations numerically along the curve 
we see that there is an excellent agreement be- 
tween LB theory and the simulations, Fig. 2. If the 
macroscopic oscillations are suppressed and the paramet- 
ric resonances are not excited, the system is able to relax 
to a quasi-ergodic equilibrium permitted by the Vlasov 
dynamics. 

To make clear the role of parametric resonances in er- 
godicity breaking, in Fig. 3a we plot the Poincare section 
of a set of non-interacting test particles, which at t = 
are distributed in accordance with Eq.Q. The motion 
of each particle is governed by Eq.Q with M x deter- 
mined by Eqs. |7]) and Q. The position and momentum 
of each particle are plotted when magnetization is at its 
minimum. We see that if the energy and the initial mag- 
netization lie on the generalized virial curve — point (B) 
of Fig. 1 — particle trajectories are completely regular. 
However, when initial conditions do not coincide with the 
generalized virial curve — point (T) of Fig. 1 — para- 
metric resonances appear and dynamics becomes chaotic. 
Particles enter in resonance with the oscillations of the 
mean-field potential, gaining sufficient energy to move 
into statistically improbable regions of the phase space. 
The Poincare section of test particle dynamics is remark- 
ably similar to the final stationary distribution obtained 
using the complete N-body molecular dynamics simula- 
tion of the HMF, Fig. 3. Eq. ^ can also be used to 
calculate the period of the first oscillation of M(t). For 
example, for point (T) of the phase diagram Fig. [TJ we 
find the period to be T = 5.0, while the full molecular 
dynamics simulations gives T = 5.4. For point (C) we 
find T — 3.85, while the simulations give T = 3.82. 

In conclusion, we have studied the mechanism respon- 
sible for the ergodicity breaking in systems with long- 
range interactions. Ergodicity breaking and the paramet- 
ric resonances are intimately connected. If the macro- 
scopic oscillations — and the resulting resonances are 
suppressed — the system is able to relax to a quasi- 
ergodic stationary state. However, when the paramet- 
ric resonances are excited, some particles are ejected to 
statistically improbable regions of the phase space, at 
the same time as the oscillations are damped out. The 
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FIG. 2. The angle and velocity distribution functions cor- 
responding to the initial conditions described by points (A), 
(B) and (C) of Fig. 1, respectively. Symbols are the results of 
molecular dynamics simulations and solid curves are the pre- 
dictions of LB theory. The simulated distribution functions 
for the point (B), lying on the generalized virial curve, are 
in excellent agreement with the predictions of the LB theory 
[panels (c) and (d)], demonstrating that the dynamics along 
the generalized virial curve is quasi-ergodic. On the other 
hand, the distribution functions for points (A) and (C) devi- 
ate significantly from the predictions of LB theory — showing 
that away from the generalized virial curve, ergodicity is bro- 
ken [panels (a), (b) and (e), (f)]. 



process of continuous particle ejection, and the result- 
ing decrease of macroscopic oscillations of the envelope, 
leads to the formation of a static mean-field potential 
and to asymptotically integrable dynamics. Once the in- 
tegrability of the equations of motion is achieved, the 
ergodicity becomes irreversibly broken. Unlike for par- 
ticles with short-range interaction potentials, ergodicity 
is the exception rather than the rule for systems with 
long-range forces — it can only be observed if the initial 
distribution function satisfies the generalized virial con- 
dition derived in this Letter. Finally we note, that since 
the stationary distribution must satisfy the virial con- 
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FIG. 3. Poincare sections of test particles and snapshots of 
the phase space obtained using molecular dynamics simula- 
tion once the system has relaxed to qSS. Panels (a) and (b) 
correspond to the initial condition lying on the generalized 
virial curve, point (B) of Fig. 1. In this case the test particle 
dynamics is completely regular, and the stationary particle 
distributions are well described by LB theory. Panels (c) and 
(d) correspond to the initial conditions slightly off the virial 
curve, point (T) of Fig. 1. Even though we have moved 
only a little from the virial curve, we see the appearance of 
resonant islands and the dynamics of some of the test parti- 
cles becoming chaotic. Such resonances drive some particles 
of the HMF to statistically improbable — from the point of 
view of the Boltzmann-Gibbs and LB statistical mechanics 
— regions of the phase space. Once the envelope oscillations 
are damped out, particle dynamics becomes completely inte- 
grable, and there is no mechanism for the resonant particles 
to equilibrate with the rest of the distribution. Thus, away 
from the generalized virial curve, ergodicity becomes broken. 



dition and the energy must be conserved, Eq.([9| allows 
us to predict the magnetization to which the system will 
evolve for initial conditions lying inside the ferromagnetic 
region. For example, point (A) of Fig. 1 which has initial 
magnetization and energy M$ = 0.74 and u = 0.55, will 
evolve to a final stationary state with M = 0.56; while 
the point (C) with Mo = 0.74 and u = 0.25, will evolve 
to a final stationary state with M = 0.86, which are pre- 
cisely the values obtained using the molecular dynamics 
simulations. 
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